practical 3

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

library(sf)
library(here)
st_layers("~/Documents/_UCL_grad school/UCL/Term 1 /GIS CASA/casa_0005_coursework/week 3/gadm36_AUS_gpkg/gadm36_AUS.gpkg")

Driver: GPKG 
Available layers:

Ausoutline <- st_read("~/Documents/_UCL_grad school/UCL/Term 1 /GIS CASA/casa_0005_coursework/week 3/gadm36_AUS_gpkg/gadm36_AUS.gpkg", 
                      layer='gadm36_AUS_0')

Reading layer `gadm36_AUS_0' from data source 
  `/Users/christine/Documents/_UCL_grad school/UCL/Term 1 /GIS CASA/casa_0005_coursework/week 3/gadm36_AUS_gpkg/gadm36_AUS.gpkg' 
  using driver `GPKG'
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 112.9211 ymin: -55.11694 xmax: 159.1092 ymax: -9.142176
Geodetic CRS:  WGS 84

print(Ausoutline)

Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 112.9211 ymin: -55.11694 xmax: 159.1092 ymax: -9.142176
Geodetic CRS:  WGS 84
  GID_0    NAME_0                           geom
1   AUS Australia MULTIPOLYGON (((158.6928 -5...

st_crs(Ausoutline)$proj4string

[1] "+proj=longlat +datum=WGS84 +no_defs"

library(raster)
library(terra)
jan<-terra::rast("~/Documents/_UCL_grad school/UCL/Term 1 /GIS CASA/casa_0005_coursework/week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_01.tif")
# have a look at the raster layer jan
jan

class       : SpatRaster 
dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : wc2.1_5m_tavg_01.tif 
name        : wc2.1_5m_tavg_01 
min value   :          -46.042 
max value   :           34.065 

#average temperature january (01 month)
plot(jan)

newproj<-"ESRI:54009"
# get the jan raster and give it the new proj4
pr1 <- jan %>%
  terra::project(., newproj)
plot(pr1)

#ok that was just testing projections

pr1 <- pr1 %>%
  terra::project(., "EPSG:4326")
plot(pr1)

# look in our folder, find the files that end with .tif and 
library(fs)
dir_info("week 3/") 

library(tidyverse)
listfiles<-dir_info("week 3/wc2.1_5m_tavg/") %>%
  filter(str_detect(path, ".tif")) %>%
  dplyr::select(path)%>%
  pull()

#have a look at the file names 
listfiles

week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_01.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_02.tif 
week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_03.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_04.tif 
week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_05.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_06.tif 
week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_07.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_08.tif 
week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_09.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_10.tif 
week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_11.tif week 3/wc2.1_5m_tavg/wc2.1_5m_tavg_12.tif 

worldclimtemp <- listfiles %>%
  terra::rast()
  
#have a look at the raster stack
worldclimtemp

class       : SpatRaster 
dimensions  : 2160, 4320, 12  (nrow, ncol, nlyr)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : wc2.1_5m_tavg_01.tif  
              wc2.1_5m_tavg_02.tif  
              wc2.1_5m_tavg_03.tif  
              ... and 9 more source(s)
names       : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ... 
min values  :     -46.042,     -44.800,     -57.986,     -64.200,     -64.829,     -64.395, ... 
max values  :      34.065,      32.908,      33.081,      34.277,      36.299,      38.458, ... 

worldclimtemp$Jan

class       : SpatRaster 
dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : wc2.1_5m_tavg_01.tif 
name        :     Jan 
min value   : -46.042 
max value   :  34.065 

#create data frame of cities in australia

site <- c("Brisbane", "Melbourne", "Perth", "Sydney", "Broome", "Darwin", "Orange", 
          "Bunbury", "Cairns", "Adelaide", "Gold Coast", "Canberra", "Newcastle", 
          "Wollongong", "Logan City" )
lon <- c(153.03, 144.96, 115.86, 151.21, 122.23, 130.84, 149.10, 115.64, 145.77, 
         138.6, 153.43, 149.13, 151.78, 150.89, 153.12)
lat <- c(-27.47, -37.91, -31.95, -33.87, 17.96, -12.46, -33.28, -33.33, -16.92, 
         -34.93, -28, -35.28, -32.93, -34.42, -27.64)
#Put all of this inforamtion into one list 
samples <- data.frame(site, lon, lat, row.names="site")
# Extract the data from the Rasterstack for all points 
AUcitytemp<- terra::extract(worldclimtemp, samples)

#create a subset for perth 
library(tidyverse)
#define where you want the breaks in the historgram
userbreak<-c(8,10,12,14,16,18,20,22,24,26)

# remove the ID and site columns
Perthtemp <- Aucitytemp2 %>%
  filter(site=="Perth")

t<-Perthtemp %>%
 dplyr::select(Jan:Dec)

  hist((as.numeric(t)), 
     breaks=userbreak, 
     col="red", 
     main="Histogram of Perth Temperature", 
     xlab="Temperature", 
     ylab="Frequency")

histinfo <- as.numeric(t) %>%
  as.numeric()%>%
  hist(.)

histinfo

$breaks
[1] 12 14 16 18 20 22 24 26

$counts
[1] 2 2 2 1 1 2 2

$density
[1] 0.08333333 0.08333333 0.08333333 0.04166667 0.04166667 0.08333333 0.08333333

$mids
[1] 13 15 17 19 21 23 25

$xname
[1] "."

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

AusoutSIMPLE <- Ausoutline %>%
  st_simplify(., dTolerance = 1000) %>%
  st_geometry()%>%
  plot()

print(Ausoutline)

Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 112.9211 ymin: -55.11694 xmax: 159.1092 ymax: -9.142176
Geodetic CRS:  WGS 84
  GID_0    NAME_0                           geom
1   AUS Australia MULTIPOLYGON (((158.6928 -5...

crs(worldclimtemp)

[1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

Austemp <- Ausoutline %>%
  # now crop our temp data to the extent
  terra::crop(worldclimtemp,.)

# plot the output
plot(Austemp)

exactAus<-terra::mask(Austemp, Ausoutline)
#subset using the known location of the raster
hist(exactAus[[3]], col="red", main ="March temperature")

exactAusdf <- exactAus %>%
  as.data.frame()

library(ggplot2)
# set up the basic histogram
gghist <- ggplot(exactAusdf, 
                 aes(x=Mar)) + 
  geom_histogram(color="black", 
                 fill="white")+
  labs(title="Ggplot2 histogram of Australian March temperatures", 
       x="Temperature", 
       y="Frequency")
# add a vertical line to the hisogram showing mean tempearture
gghist + geom_vline(aes(xintercept=mean(Mar, 
                                        na.rm=TRUE)),
            color="blue", 
            linetype="dashed", 
            size=1)+
  theme(plot.title = element_text(hjust = 0.5))

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

squishdata<-exactAusdf%>%
  pivot_longer(
  cols = 1:12,
  names_to = "Month",
  values_to = "Temp"
)

twomonths <- squishdata %>%
  # | = OR
  filter(., Month=="Jan" | Month=="Jun")

meantwomonths <- twomonths %>%
  group_by(Month) %>%
  summarise(mean=mean(Temp, na.rm=TRUE))

meantwomonths

ggplot(twomonths, aes(x=Temp, color=Month, fill=Month)) +
  geom_histogram(position="identity", alpha=0.5)+
  geom_vline(data=meantwomonths, 
             aes(xintercept=mean, 
                 color=Month),
             linetype="dashed")+
  labs(title="Ggplot2 histogram of Australian Jan and Jun
       temperatures",
       x="Temperature",
       y="Frequency")+
  theme_classic()+
  theme(plot.title = element_text(hjust = 0.5))

data_complete_cases <- squishdata %>%
  drop_na()%>% 
  mutate(Month = factor(Month, levels = c("Jan","Feb","Mar",
                                          "Apr","May","Jun",
                                          "Jul","Aug","Sep",
                                          "Oct","Nov","Dec")))

# Plot faceted histogram
ggplot(data_complete_cases, aes(x=Temp, na.rm=TRUE))+
  geom_histogram(color="black", binwidth = 5)+
  labs(title="Ggplot2 faceted histogram of Australian temperatures", 
       x="Temperature",
       y="Frequency")+
  facet_grid(Month ~ .)+
  theme(plot.title = element_text(hjust = 0.5))

library(plotly)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:raster’:

    select

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

# split the data for plotly based on month

jan <- squishdata %>%
  drop_na() %>%
  filter(., Month=="Jan")

jun <- squishdata %>%
  drop_na() %>%
  filter(., Month=="Jun")

# give axis titles
x <- list (title = "Temperature")
y <- list (title = "Frequency")

# set the bin width
xbinsno<-list(start=0, end=40, size = 2.5)

# plot the histogram calling all the variables we just set
ihist<-plot_ly(alpha = 0.6) %>%
        add_histogram(x = jan$Temp,
        xbins=xbinsno, name="January") %>%
        add_histogram(x = jun$Temp,
        xbins=xbinsno, name="June") %>% 
        layout(barmode = "overlay", xaxis=x, yaxis=y)

ihist

# mean per month
meanofall <- squishdata %>%
  group_by(Month) %>%
  summarise(mean = mean(Temp, na.rm=TRUE))

# print the top 1
head(meanofall, n=1)

# standard deviation per month
sdofall <- squishdata %>%
  group_by(Month) %>%
  summarize(sd = sd(Temp, na.rm=TRUE))

# maximum per month
maxofall <- squishdata %>%
  group_by(Month) %>%
  summarize(max = max(Temp, na.rm=TRUE))

# minimum per month
minofall <- squishdata %>%
  group_by(Month) %>%
  summarize(min = min(Temp, na.rm=TRUE))

# Interquartlie range per month
IQRofall <- squishdata %>%
  group_by(Month) %>%
  summarize(IQR = IQR(Temp, na.rm=TRUE))

# perhaps you want to store multiple outputs in one list..
lotsofstats <- squishdata %>%
  group_by(Month) %>%
  summarize(IQR = IQR(Temp, na.rm=TRUE), 
            max=max(Temp, na.rm=T))

# or you want to know the mean (or some other stat) 
#for the whole year as opposed to each month...

meanwholeyear=squishdata %>%
  summarize(meanyear = mean(Temp, na.rm=TRUE))